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Abstract. We calculate energy spectra and mass composition of cosmic rays accelerated by the galactic pop- 
ulation of pulsars during their radio and gamma-ray phase. It is assumed that a significant part of the pulsar 
rotational energy is lost on acceleration of iron nuclei extracted from the surface of the neutron star. The nuclei 
are accelerated, at first when passing the outer gap of the inner pulsar magnetosphere and later in the pulsar wind 
zone, to energies corresponding to 50% of the total potential drop through the polar cap region of the neutron 
star. We calculate energy spectra of the nuclei injected from the pulsar wind nebulae into the Galaxy including 
different energy loss processes of nuclei during their propagation in the pulsar magnetosphere and the expanding 
nebula, their fragmentation in collisions with radiation and matter, adiabatic energy losses, and escape from the 
nebula. Several models proposed for the distribution of the initial parameters of the galactic pulsar population are 
considered. It is shown that the best description of the observed cosmic ray spectrum and the mass composition 
between a few 10 15 eV and a few 10 18 eV is obtained for the model B of Lorimer et al., in which the logs of 
initial pulsar periods and surface magnetic fields are given by the Gaussian distributions with the average values 
of < logP[ms] >= 2.6 and < logB[G] >= 12.3, respectively. In order to supply the cosmic rays into the Galaxy 
with the required rate, the product of the efficiency of cosmic ray acceleration by pulsars and their birth rate 
should be about 10~ 2 yr -1 . 

Key words, supernova remnants: pulsars: general - cosmic rays 



1. Introduction 

The origin of cosmic rays (CRs) with energies around and 
above ~ 10 15 eV, the knee region, is still far from be- 
ing understood. It is widely accepted that the lower en- 
ergy CRs, below the knee, are accelerated at large scale 
shock waves created in the interstellar medium by super- 
nova explosions. It is argued that in some favorite condi- 
tions the particles might be also accelerated above 10 15 
eV by this same supernova shock mechanism if, e.g. mul- 
tiple supernovae exploded in the star forming regions - 
superbubbles (Cesarsky & Montmerle 1983, Bykov 2001), 
a supernova explodes in its dense wind (Biermann 1993, 
Stanev et al. 1993), the magnetic field at the shock front is 
non-linearly amplified by CRs (Lucek & Bell 2000, Bell & 
Lucek 2001, Drury et al. 2003), or superheavy nuclei are 
accelerated (Horandel 2003). The acceleration of particles 
to such energies might also occur in a large scale terminal 
shock of the galactic wind (Jokipii & Morrill 1985). Other 
models propose that the observed cosmic rays has extra- 
galactic origin. For example, Protheroe & Szabo (1992) 
consider their acceleration in the shocks in the accre- 
tion flows in active galactic nuclei, and Vietri (1995) 
and Waxman (1995) argue for acceleration on relativis- 
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tic shocks of gamma-ray bursts (GRBs). This last idea 
has been considered recently in a more detail by Wick 
et al. (2004), who propose the origin of cosmic rays from 
~ 10 14 eV up to the highest energies in the GRBs includ- 
ing the events which occurred in our own Galaxy with 
typical rate of (3 — 10) x 10 -6 per yr. For the extensive 
review of recent models of the cosmic ray acceleration we 
refer to Horandel (2004). 

Another acceleration mechanism was proposed long 
time ago, soon after the pulsar discovery. It postulated ac- 
celeration of CRs above the knee region in the pulsar wind 
zone in large amplitude electromagnetic waves generated 
by rotating neutron star (Ostriker & Gunn 1969, Karakula 
et al. 1974). The contribution of particles accelerated by 
pulsars to the observed cosmic ray spectrum has been 
more recently discussed in detail by Cheng & Chi (1996), 
Bednarek & Protheroe (2002), Giller & Lipski (2002). For 
example, Bednarek & Protheroe (2002) estimate the con- 
tribution of heavy nuclei accelerated in the pulsar outer 
gaps (Cheng et al. 1986) to energies above the knee region. 
Giller & Lipski (2002) derive the initial parameters of the 
pulsar population inside the Galaxy required to explain 
the observed shape of the CR spectrum and its inten- 
sity. Note that also CRs with extremely high energies, i.e. 
above the ankle at ~ 10 18 eV, might be accelerated in the 
wind regions of pulsars with super strong surface mag- 
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netic fields, objects called magnetars (Blasi et al. 2000, 
Arons 2003). 

The models postulating acceleration of CRs between 
the knee and the ankle locally in our Galaxy obtained 
strong support by recent observation of the cosmic ray 
anisotropies from the direction of the Galactic Centre 
and the Cygnus region (Hayashida et al. 1999, Bellido et 
al. 2001). The AGASA experiment (Hayashida et al. 1999) 
reported the excesses of CRs in the narrow energy range 
close to ~ 10 18 eV with the significance of 4.5<r and 3.9a 
for the Galactic Centre and the Cygnus region, respec- 
tively. 

Another information on the origin of CRs can be 
reached by measuring its mass composition. Most of the 
experiments suggest that the average mass composition 
of CRs is relatively light below the knee (e.g. Burnett et 
al. 1990, Apanasenko et al. 2001) and becomes heavier 
above the knee (e.g. Glasmacher et al. 1999, Arqueros et 
al. 2000, Fowler et al. 2001, Shirasaki et al. 2001, Ave et 
al. 2003). Recent preliminary results obtained from the 
analysis of the KASCADE data confirm heavier compo- 
sition above the knee and show evidences of the rigidity 
dependent break in the spectra of specific groups of nu- 
clei (e.g. Roth et al. 2003). The mass composition seems 
to decline again to lighter between 10 17 — 10 18 eV, i.e. 
at the region of the second knee, according to some ex- 
periments (Bird et al. 1993, Yoshida et al. 1995, Dawson 
et al. 1998), or may be still quite heavy in this energy 
range as it is obtained in the analysis of the Haverah Park 
data which show bi-modal composition composed of one 
third of protons and two thirds of iron (Abu-Zayyad et 
al. 2001). Recent analysis of the Vulcano Range data also 
suggest heavy composition composed of 88% of the iron 
nuclei in the bi-modal proton-iron composition at a me- 
dian energy of 10 18 eV (Dova et al. 2003). The change of 
the cosmic ray mass composition at the knee region is of- 
ten explained as a result of rigidity dependent acceleration 
causing heavy nuclei start to dominate at higher energies 
(Peters 1961). The heavy nuclei enrichment of the galac- 
tic CRs might be also due to propagation effects in the 
galactic disk and halo (e.g. Maurin et al. 2003). Moreover, 
a contribution of another type of sources to CRs above 
the knee region, e.g. pulsars, may also increase its average 
mass. 

In this paper we consider the last possibility in more 
details. We calculate the energy spectra of different types 
of nuclei injected by the pulsars into the Galaxy, apply- 
ing a model for their injection and propagation in the ex- 
panding pulsar wind nebula (PWNa) surrounding a young 
pulsar. The model takes into account the energy losses 
and escape conditions of nuclei during the expansion of 
the nebula. Similar model has been recently applied for 
modelling the 7-ray emission from the PWNe (Bednarek 
& Bartosik 2003). We estimate the contribution of nu- 
clei, escaping from the PWNe, to the observed cosmic ray 
spectrum above the knee for different models of the pul- 
sar population in the Galaxy. The results of calculations 



are compared with the reports on the mass composition 
in this energy region. 

2. The energy spectra of nuclei injected from the 
pulsar wind nebulae 

In this section we calculate spectra of nuclei accelerated 
by young pulsars and injected from the pulsar wind nebu- 
lae into the Galactic medium. We define the acceleration 
mechanism of nuclei, consider the effects of their propaga- 
tion in the expanding nebula, and discuss the results for 
different models of the galactic pulsar population. 

2.1. Acceleration of nuclei by the pulsar 

It is likely that rotating magnetospheres of neutron stars 
can accelerate not only leptons but also heavy nuclei, ex- 
tracted from positively charged polar cap regions. In fact 
different aspects of high energy phenomena around pul- 
sars may need the presence of heavy nuclei. For exam- 
ple, a flow of iron nuclei, coexisting with an outflowing 
electron-positron plasma, is postulated in order to explain 
the change in the drift direction of the radio subpulses 
(Gil et al. 2003). The presence of heavy nuclei can also 
explain morphological features of the Crab Nebula and 
the appearance of extremely energetic leptons inside the 
nebula, accelerated as a result of resonant scattering of 
positrons and electrons by heavy nuclei (Hoshino et al. 
1992, Gallant & Arons 1994). From normalization to the 
observations of the Crab pulsar, Arons and collaborators 
(e.g. see Arons 1998) postulate that the Lorentz factors 
of iron nuclei, accelerated somewhere in the inner magne- 
tosphere and the pulsar wind zone and, injected into the 
pulsar wind nebula should be, 

7 jFe « xZe<P opcn /m iC 2 « 8 x l0 9 X B 12 P m 2 , (1) 

where m; and Ze are the mass and charge of the iron nu- 
clei, c is the velocity of light, and $ pcn = y Lmt / c is 
the total electric potential drop across the open magne- 
tosphere, L m t is the rate of rotational energy lost by the 
pulsar, B = \Q 12 B\2G is the surface magnetic field of the 
pulsar, and P — 10 _3 P ms is the pulsar period. Due to the 
rotational energy losses on emission of dipole radiation, 
the pulsar period evolves in time according to 

^ s = Plms + 2 x 10- 9 tB 2 12 , (2) 

where i"b,ms is the initial period of the pulsar, and t is in 
seconds. Arons and collaborators argue that the acceler- 
ation factor x is n °t far from unity, x = 0-5 is taken in 
the following calculations. Moreover, we assume that these 
nuclei take significant part, £, of the pulsar rotational en- 
ergy. The unknown value of £ times the pulsar birth rate 
in the Galaxy, 77, can be obtained from the comparison 
of contribution of particles injected by pulsars to the ob- 
served cosmic ray flux. Eq. (JTJ postulates that the pul- 
sar at a specific time accelerates nuclei monoenergetically. 
However the energies of freshly produced nuclei change in 
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time due to the change of the pulsar period caused by its 
rotational energy losses. 

We assume that a pulsar can extract iron nuclei from 
its surface and accelerate them to high energies during the 
radio phase. Then the products of efficient leptonic cas- 
cades, occurring in the inner pulsar magnetosphere, im- 
pinge on the polar cap region allowing extraction of the 
iron nuclei from the surface. 

The iron nuclei extracted from the surface of the pul- 
sar are accelerated in the outer gaps of the inner pul- 
sar magnetosphere (Cheng et al. 1986), suffering partial 
photo-disintegrations in collisions with nonthermal radi- 
ation produced in leptonic cascades. The details of the 
acceleration model of nuclei in the outer gaps are given in 
section 3.1 of Bednarek & Protheroe (2002). The level of 
photo-disintegration of nuclei depends on the pulsar pa- 
rameters. To calculate the number of nucleons extracted 
from nuclei during propagation through the outer gap 
we determine the outer gap radiation field (Bednarek & 
Protheroe 2002), by interpolation/extrapolation of the 
radiation fields in the Crab type pulsars, calculated by 
Ho (1989), and for the Vela type pulsars, based on the 
observed change of the 7-ray luminosity (Ruderman & 
Cheng 1988). Nuclei and protons from their disintegra- 
tion are additionally accelerated in the pulsar wind zone 
to energies postulated by Arons (1998), see Eq. 1. They 
are injected into the pulsar wind nebula. Neutrons from 
the disintegration of nuclei move ballistically through the 
nebula decaying inside or outside it depending on their 
Lorentz factor. At the early stage of expansion, when the 
nebula is dense enough, nuclei suffer further disintegra- 
tions in collisions with the matter (Bednarek & Protheroe 
2002). These processes are discussed in the next subsec- 
tion. 

2.2. Escape of nuclei from expanding nebula 

The nuclei, injected by the pulsar, propagate in the ex- 
panding nebula with parameters changing drastically in 
time. In order to take properly into account different ef- 
fects on their propagation (collisions with matter, diffusion 
inside the nebula, escape from the nebula) , we have to as- 
sume a time dependent model for the expanding nebula 
taking into account not only the initial parameters of the 
expanding envelope but also the energy supplied by the 
pulsar. The model has to determine such basic parameters 
as: the expansion velocity of the nebula, the pulsar wind 
shock radius, the outer radius of the nebula, the average 
density of matter and magnetic field strength inside the 
nebula, and others. The details of such a simple model for 
the PWNe are described in our previous paper (Bednarek 
& Bartosik 2003). We repeat here its main features. 

The evolution of a supernova remnant, containing an 
energetic pulsar, is described according to the picture con- 
sidered by Ostriker & Gunn (1971) and Rees & Gunn 
(1974). Let us denote the initial expansion velocity of the 
bulk matter in supernova envelope by Vq,sn and its ini- 



tial mass by Mo,sn- The expansion velocity can increase 
due to the additional supply of energy to the nebula by 
the pulsar. It can also decrease due to the accumulation 
of the surrounding matter. We take these processes into 
account to determine the radius of the nebula at a specific 
time, t, by using the energy conservation, 



M SN (*)V s 2 N (t) _ M ^V 2 

SN 



+ L P ui- 



ncb: 



(3) 



where £ pu i-ncb is the part of the energy supplied by the 
pulsar to the nebula in the form of magnetic field and the 
energy transfered from relativistic hadrons to the nebula 
due to adiabatic energy losses. The rest of the rotational 
energy lost by the pulsar is almost immediately radiated 
by accelerated leptons. The energy of relativistic hadrons 
which is injected into the interstellar medium from the 
nebula is equal to the difference between the energy gained 
by hadrons during acceleration process, assumed equal to 
£L ro t, and the energy lost by them to the nebula during its 
adiabatic expansion. The part of energy lost by relativis- 
tic hadrons during adiabatic expansion of the nebula has 
been calculated numerically since the process of diffusion 
of hadrons inside the nebula has to be taken into account. 

The mass of the expanding nebula increases (by sweep- 
ing the surrounding medium) according to the formula 



M SN (t) = M , SN + -np SUI R^ eh (t), 



(4) 



where p sm is the density of the surrounding medium and 
i?Nob is the outer radius of the expanding envelope at the 
time, t, which depends on the expansion history of the 
nebula, 

i?Ncb = / V SN (t')dt'. (5) 
Jo 

The expansion velocity of the nebula, Vsn(£), and the av- 
erage density of matter inside it, 



PNcb = 3M SN (t)/4irR^ eh (t), 



(6) 



at an arbitrary time have been found by solving numeri- 
cally the above set of Eqs. (3-6). Note however, that the 
average density of matter, which is very useful for our sim- 
ple model, may not correspond to the real distribution of 
matter inside the nebula specially at the later phase of ex- 
pansion when the part of matter can form dense filaments. 
For example the matter inside the Crab Nebula is accumu- 
lated in filaments containing 4.6 ±1.8 Mq of ionized and 
neutral gas (Fesen et al. 1997) with density <~ 500 cm~ 3 
(Davidson & Fesen 1985). No significant gas has been con- 
firmed around the Crab Nebula (Fesen et al. 1997). Atoyan 
& Aharonian (1996) suggest that relativistic hadrons can 
be efficiently trapped by these dense filaments. If this is 
the case then the efficiency of interactions of relativistic 
nuclei with the matter can be enhanced in respect to the 
interactions with the uniformly distributed matter inside 
the nebula. 

The pulsar loses energy in the form of a relativistic 
wind extending up to the shock at a distance R s h- At 
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Fig. 1. Energy spectra of different types of nuclei escaping from a supernova nebula into the galactic medium. The 
nuclei are accelerated by the pulsar with the surface magnetic field B = 4x 10 12 G and different initial periods (a) 1 
ms, (b) 10 ms, and (c) 100 ms. The pulsar has been born during the supernova explosion with the envelope mass of 
AMq expanding with the initial velocity 2000 km s _1 . The total spectrum of nuclei is marked by the thick full curve. 
Specific curves show the spectra of nuclei with the mass numbers A = 56 (dot-dot-dot-dashed), 41-55 (dotted), 11-40 
(dot-dashed), 2-10 (dashed), protons and protons from decay of neutrons (thin full). The spectra of nuclei for pulsars 
with different initial periods 1 ms (full curve), 10 ms (dot-dashed), 10 2 ms (dotted), and 10 3 ms (dot-dot-dot-dashed) 
are shown in (d). The case of the pulsar with 1 ms initial period and the expansion velocity of 10 km s _1 is shown 
by the dashed curve. 



this distance, the pressure of the wind is balanced by the 
pressure of the expanding nebula. Rees & Gunn (1974) 
estimate the location of this shock as a function of time 
by comparing the pulsar wind energy flux, determined by 
L mt (see Eq. JHJ, with the pressure of the outer nebula, 
determined by the supply of energy to the nebula by the 
pulsar over the whole of its lifetime, 

-^rot(^) -^pul — neb /— \ 

47r n2 ~ 4 n3 ' ' > 

^ ntX sh C 3 ''"Neb 

and 

L Iot (t) = B 2 s Rtn A /&c z « 2.5 x 10 43 B 2 2 P ms 4 erg s" 1 , (8) 

where R s and B s are the radius of the pulsar and its sur- 
face magnetic field, £1 = 2ir/P. The pulsar period changes 
according to Eq. It is assumed in this formula that 

i?sh "C -RNeb- 

Knowing how the magnetic field depends on the dis- 
tance from the pulsar in the pulsar wind zone, we can 



estimate the strength of the magnetic field at the shock 
region from 

where a is the ratio of the magnetic energy flux to the 
particle energy flux from the pulsar at the location of the 
pulsar wind shock. There is observational evidence that a 
evolves in time. It is much less than one for very young 
nebulae (of the Crab type), and increases with age be- 
ing of the order of one for nebulae of the Vela type. The 
value of a for a specific pulsar is determined by the effi- 
ciency of conversion of magnetic energy into the particles 
energy in the pulsar wind zone. The evolution of a with 
the parameters of the pulsar is estimated by interpolating 
between the values for the Crab pulsar, equal to ~ 0.003 
(Kennel & Coroniti 1984), and for the Vela pulsar, equal 
to ~ 1 (Hclfand ct al. 2001). We assume that the mag- 
netic field strength in the volume of the nebula depends 



W. Bednarek & M. Bartosik: Cosmic rays from Galactic pulsars 



on the distance from its center as described by the formu- 
lae derived in the magnetohydrodynamic model for the 
PWNe considered by Kennel & Coroniti (1984). The ge- 
ometry of the magnetic field inside the nebula may have 
strong effect on the diffusion of nuclei inside it. Recent 
3-dimensional modelling of the Chandra observations of 
the Crab Nebula suggests that the regular component of 
the magnetic field should be comparable to the turbu- 
lent component (Shibata et al. 2003). In our calculations 
the diffusion in turbulent component is only considered. 
The presence of regular toroidal magnetic field should slow 
down the diffusion of nuclei through the nebula. 

It is assumed that nuclei injected into the nebula at a 
specific time ii„j, escape from it at the time i esc , if their dif- 
fusion distance in the turbulent magnetic field of the neb- 
ula, -Rdiff: is equal to the dimension of the nebula, R^ch, 
at the time t esc . The diffusion distance in the magnetic 
field of the nebula is obtained by integration (Bednarek & 
Protheroe 2002), 




where the diffusion coefficient is taken to be D = R^c/3, 
and i?L is the Larmor radius of nuclei depending on the 
distance from the center of the nebula. i?diff and t csc have 
been calculated numerically. We include also the adiabatic 
losses of nuclei during the period of their propagation in- 
side the nebula. Due to this losses and interaction of nuclei 
with the matter the energy of a nucleus at a specific time 
is 

where E(ti n j) is its energy at the time ti n j, K is the inelas- 
ticity coefficient for collisions of nuclei with the matter, 
and r is the optical depth calculated based on the known 
density of matter inside the nebula (see Eq. 7). The de- 
tails of these calculations are given in sections 4.2 and 4.3 
in Bednarek & Protheroe (2002). 

Considering all the processes discussed above, we cal- 
culate the energy spectra of different groups of nuclei es- 
caping from the nebula into the Galaxy for arbitrary ini- 
tial parameters of the pulsars. Protons from decay of neu- 
trons, released from heavy nuclei during their propagation 
through the radiation field of the outer gap of the pulsar 
and later as a result of collisions with the matter inside the 
nebula, are also included in these calculations. Neutrons 
extracted from nuclei in collisions with the outer gap ra- 
diation field have relatively low Lorentz factors since their 
parent nuclei are pre-accelerated in the outer gap to en- 
ergies lower than that achieved later in the pulsar wind 
zone. On the other hand, neutrons extracted from nuclei 
in collisions with the matter have typically higher Lorentz 
factors than protons since, moving balistically through the 
PWNa, do not lose energy on the adiabatic expansion. 
We assume that neutrons decaying within 3 kpc from the 
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galactic disk are captured by the magnetic field of the 
galactic halo and contribute to the galactic CRs. As an 
example, Figs. 1 show the spectra of different types of 
nuclei injected into the Galaxy for the case of a super- 
nova with initial parameters derived for the Crab Nebula, 
i.e. the mass of the supernova envelope equal to 4 M Q , 
the initial expansion velocity 2000 km s _1 , and the den- 
sity of surrounding medium equal to 0.3 cm -3 . Figs, la-c 
show spectra of the escaping nuclei within selected range 
of mass numbers for pulsars with initial periods 1 ms, 
10 ms, and 100 ms, and a typical surface magnetic field 
of 4 x 10 12 G. It is clear that nuclei injected by pulsars 
with shorter initial periods suffer more severe fragmenta- 
tions. The total spectra of these secondary nuclei can be 
even dominated in some range of energies by secondary 
protons extracted from heavy nuclei and from decay of 
neutrons (thin full curves). The total spectra of the nu- 
clei above ~ 10^ 7-8 ' GeV have the spectral index not far 
from -2 extending through about two decades in energy. 
The low energy bump in the spectra below <~ 10 6 GeV 
is due to the neutrons extracted from nuclei during their 
fragmentation in the radiation field of the outer gap of 
the pulsar. Fig. Id compares the total spectra of nuclei in- 
jected by pulsars with different initial periods and initial 
expansion velocities of the supernova envelopes. Pulsars 
with longer initial periods should be able to contribute to 
lower energy range in the total spectrum of nuclei. On the 
other hand, larger initial expansion velocity of the enve- 
lope allows the nuclei to escape at earlier times because 
the magnetic field inside the envelope drops more rapidly. 
Therefore, these nuclei escape with larger energies due to 
their smaller energy losses on the adiabatic expansion. 

2.3. Pulsar population within the Galaxy 

The spectrum of nuclei injected inside the Galaxy can be 
obtained by summing up over all population of pulsars 
inside the Galaxy. However, the initial parameters of the 
pulsars are not precisely known. Therefore, we consider 
a few different models proposed in the literature. They 
differ in distributions of the surface magnetic fields and 
the initial periods of the new born pulsars. The following 
models are: 

1. The surface magnetic fields of the pulsars are described 
by the distribution derived by Narayan (1987), 

dN/d(logB) w 0.065/(B[10 12 G]) year" 1 , (12) 

for B > 2 x 10 12 G. The distribution of pulsars with 
the surface magnetic fields below 2 x 10 12 G bases on 
Fig. 13 in Narayan (1987). All pulsars are born with 
the fixed initial period Po = 40 ms (van der Swaluw & 
Wu 2001). 

2. The surface magnetic fields of the pulsars as in model 
(i) . The initial pulsar periods are correlated with their 
surface magnetic fields, 

P [ms] = 63.7/(B[10 12 G]), (13) 
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Fig. 2. The spectra of nuclei, injected by the galactic population of pulsars into the Galaxy, calculated for different 
models of initial parameters of the pulsar population, obtained in the work by Narayan (1987) and van der Swaluw 
& Wu (2001) (model i), Narayan (1987) and Xu et al. (2001) (model ii), Model A of Lorimer ct al. (1993) (our model 
hi), and Model B of Lorimer et al. (1993) (model iv) (see text for details). The spectra of different types of nuclei are 
marked as in Figs. 1. 



as postulated by Xu et al. (2001). We apply this for- 
mula for the pulsar initial periods above 2 ms. 

3. Model A of Lorimer ct al. (1993) postulating the 
Gaussian distribution of log B of the pulsars with pa- 
rameters < logB[G] >= 12.46 and cr logB = 0.31. We 
apply this model for the range of B = io 11 - 5-13 - 3 
G. All pulsars are born with a single initial period, 
logPo[ms] — 1.35. 

4. Model B of Lorimer et al. (1993) postulating the 
Gaussian distributions of log B, with < logB[G] >= 
12.3 and <7i og B = 0.35, and for log Po, with < 
logP [ms] >= 2.6 and <J\ og B = 0.35 above 2 ms. 

For these four models of the pulsar population, we cal- 
culate the injection spectra of nuclei into the Galaxy ap- 
plying the values of the normalization coefficients, £ • 7], 
which give us the flux of calculated nuclei on the level 
comparable to that observed in the cosmic ray spectrum. 
It is assumed that the mass of nebula surrounding the pul- 
sar is equal to 3 M and expends with the velocity of 2000 
km s _1 . Model (i), calculated for £ • i] = 1/700 yr _1 , gives 
a power law spectrum of nuclei over about two decades 
above ~ 10 s GeV and very sharp cut-off at lower ener- 
gies (see Fig. 2i). The spectrum is dominated by heavy 



nuclei from the iron group. Model (ii), calculated also for 
£ • f] = 1/700 yr _1 , gives very flat spectrum at high en- 
ergies with the spectral index close to 2 (see Fig. 2ii). 
The primary iron nuclei disintegrate significantly in this 
model. Protons from decay of neutrons can contribute to 
the total spectrum at around 10 18 eV. These features are 
due to a relatively large number of pulsars with short pe- 
riods predicted by this model. Model (iii), calculated for 
£ • r/ = 1/120 yr _1 , gives the total spectrum with a rather 
narrow, strong peak at <~ 10 8 — 10 9 GeV (see Fig. 2iii). The 
total spectrum is dominated by heavy nuclei from the iron 
group. Such shape of the spectrum is due to the assump- 
tion about the fixed period of all new born pulsars which 
is probably not very realistic. Model (iv), £ • r\ = 1/120 
yr _1 , predicts the power law spectrum of nuclei between 
a few 10 7 GeV and 10 10 GeV with the index close to 3 
(see Fig. 2iv). The spectrum is dominated by heavy nuclei 
in all this energy range. In the next section we compare 
the spectra of nuclei calculated in terms of these models 
with the observed spectrum and the mass composition of 
CRs. 
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Fig. 3. The comparison of the spectrum of particles ac- 
celerated by pulsars with the parameters described by 
the model (i) (dot-dot-dot-dashed curve in (a)), model 
(ii) (dashed curve in (a)), model (iii) (dot-dot-dot-dashcd 
curve in (b)), and model (iv) (dashed curve in (b)), 
with the observed cosmic ray spectrum (thin full curve). 
It is assumed that at low energies the supernova rem- 
nants (SNR) accelerate CRs with the spectrum dN/dE oc 
E~ 2 - 7 exp(-E/10 7 GeV) (dot-dashed curve). At extremely 
high energies (EG), the spectrum has the form dN/dE oc 
E~ 2 - 7 (dotted line). 



3. Contribution to the cosmic rays in the Galaxy 

Having in hand the injection spectrum of nuclei by pulsars, 
we can estimate the flux of these nuclei inside the Galaxy 
adopting the leaky box model from the relation 



dN 



dN; 



C T e! 



dEdSdtdQ dEdt 4tt V 



gal 



(14) 



m /dEdt is the total injection spectrum of CRs 



where dN\ 

from all the pulsars in the Galaxy (calculated in Sect. 2), 
V ga i = 10 68 cm 3 is the volume of the Galaxy (a disk with 
a radius of 15 kpc and a half-thickness of 3 kpc), and 
t csc is the escape time of hadrons from the Galaxy. The 
dependence of the escape time of nuclei on their energy 
and charge is approximated by 



7"csc — 



yr 



(E GeY /Z) 



0.3 - 



(15) 



where Z is the charge of nuclei and Eq c y its energy in 
GeV. This formula is normalized to the estimated lifetime 
of protons with energies 10 18 eV equal to 4 x 10 4 years 
(see Fig. 4.18 in Berezinsky et al. 1990). 

The calculations have been performed for all four mod- 
els of the pulsar population, applying the normalization 
parameter, £-77, as reported above. Our aim is to find out if 
any of the proposed so different pulsar population models 
is able to describe general features of the observed cosmic 
ray spectrum. The comparison with the CR and calculated 
spectra are shown in Figs. 3. When fitting the observed CR 
spectrum in the broad energy range, we assume the pop- 
ular hypothesis suggesting that the shape of the observed 
cosmic ray spectrum is a combination of four different 
components. Below the knee the power law spectrum has 
the form oc E~ 2 - 7 . It is widely believed that the main con- 
tribution to this part of the spectrum comes from the par- 
ticles accelerated in the supernova shock waves resulting in 
a power law spectrum and the exponential cut-off which 
we put at 10 7 GeV, dN/dE oc E- 2 - 7 exp(-E/10 7 GeV). 
Between the knee and the second knee at a few 10 17 eV, 
the spectrum has the form E~ 3 , and between the sec- 
ond knee and the ankle oc E~ 3 - 2 . The extremely high en- 
ergy part of the cosmic ray spectrum is due to another 
component approximated by a simple power law with the 
spectral index -2.7. Particles, accelerated by the pulsars, 
dominate the spectrum between the knee and the ankle. 

However, only model (iv) gives good general consis- 
tency of the observed shape of the spectrum with the cal- 
culated one. Model (i) predicts significant deficit of parti- 
cles at ~ 10 16 eV, which has not been reported by experi- 
mentalists. This deficit might be suppressed if the pulsars 
are able to accelerate nuclei for longer time than the radio 
activity period expected from the modeling of magneto- 
spheric high energy cascade processes. Models (ii) and (iii) 
predict too flat spectra between the knee and the ankle 
and thus should be rejected. 

We also calculate the expected mass composition of 
the CRs, < In A >, in terms of the model (iv) and com- 
pare it with the measurements of the mass composition 
of CRs by different experiments (see Fig. 4). The mass 
composition of the galactic supernova component is taken 
as reported by the direct measurements at low energies 
(JACEE, RUN JOB). The same composition is also ap- 
plied to the extremely high energy component. We cal- 
culate the average mean < InA > for our model (iv) as- 
suming two models for extend of the extragalactic (EG) 
component. In the first one, EG component extends down 
to at least 3 x 10 7 GeV. In the second one, it breaks at 10 9 
GeV (see Fig. 4). The Model (iv) with EG component ex- 
tending to low energies describes well the general tendency 
reported by most experiments, i.e. the rise of < InA > 
above the knee, and gradual decrease of < InA > above 
~3x 10 17 eV. Therefore, we conclude that the model of 
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Fig. 4. The comparison of the mass composition of the 
cosmic rays, < In A >, obtained in terms of the model (iv) 
for the case of the light extragalactic component extending 
to low energies (thick dashed curve) and breaking at 10 18 
cV (thin dashed curve), with the measurements reported 
by different experiments (JACEE - Shibata et al. (1999), 
RUNJOB - Apanasenko et al. (2001), CASA-MIA - 
Glasmacher et al. (1999), KASCADE (e/m) - Ulrich et 
al. (2001), KASCADE (p/m) - Engler et al. (1999), EAS- 
TOP - Alessandro et al. (2001), FLY'S-EYE - Dawson 
et al. (1998), Akeno Al - Hayashida et al. (1995), Akeno 
A100 - Hayashida et al. (1997), Havcrah Park - Ave et 
al. (2003), Volcano Ranch - Dova et al. (2003)). The re- 
gion between the thin full curves shows < InA > obtained 
from many experiments and derived by Horandel (2004). 
< InA > equal to 1.5 is taken for the SNR contribution to 
cosmic ray spectrum as measured at low energies by di- 
rect experiments. Similar composition is also applied for 
the extragalactic (EG) component of CRs at the highest 
energies. 

the pulsar contribution to the cosmic ray spectrum can 
explain the change of the average mass above the knee, 
which is the consequence of a sudden switching off the 
pulsar acceleration mechanism of heavy nuclei due to the 
lack of their extraction from the neutron star surface. In 
our model this is due to the switching off the process of 
electromagnetic cascading in the inner pulsar magneto- 
sphere whose products heat up the polar regions of the 
neutron star. If the light extragalactic component does 
not extend to low energies but breaks at ~ 10 18 eV, then 
our model predicts rather heavy composition up to ~ 10 18 
eV as suggested by some groups (Dawson et al.1998, Ave 
et al. 2003, Dova et al. 2003). 

4. Discussion and Conclusions 

The best fit to the observed cosmic ray spectrum between 
the knee and the ankle is obtained for the model B of 
Lorimer et al. (1993) (our model iv) which postulates that 
the observed radio pulsars are born with relatively long 



initial periods, with the average value of <~ 400 ms, with 
the Gaussian distribution in log scale, and typical surface 
magnetic fields, with the average value of 2 x 10 12 G and 
also the Gaussian distribution in log scale. This average 
value of the pulsar periods at birth is similar to that one 
obtained in the analytical work by Giller & Lipski (2002), 
who got the best description of the cosmic ray spectrum 
for the gamma distribution of the initial pulsar periods 
(which is oc Pq~ for small periods) with average value of 
Pq = 500 ms and s = 3.86, and the Gaussian distribution 
for log B with average value between 10 12 — 10 13 G and 
ciogB — 0-4- However, the shape of the pulsar distribution 
for initial periods is very different from that one derived 
by Lorimer et al. (1993), and applied in our work. This is 
not surprising since, in contrast to Giller & Lipski (2002), 
we take into account the propagation and escape condi- 
tions, and adiabatic and collisional energy losses of nuclei 
during their passage throughout the pulsar wind nebula. 
These processes have strong effects on the injection spec- 
trum of nuclei into the Galaxy. From the normalization 
of the calculated spectrum of nuclei in terms of the model 
(iv) to the observed cosmic ray spectrum, we obtained the 
efficiency of conversion of the rotational energy of the pul- 
sar into the energy of iron nuclei multiplied by the pulsar 
birth rate equal to ~ 1/120 yr _1 . This value is consistent 
with the estimated pulsar birth rate <~ 1/(20 — 250 yr) 
(Lync et al. 1985; Lorimer 2003, Vranesevic et al. 2003). 
Therefore, we conclude that the efficiency of acceleration 
of the iron nuclei should be in the range ~ (0.2 — 1.0). 

The low energy break in the spectrum of nuclei injected 
by the PWNe is explained in our model by switching off 
the mechanism of extraction of the iron nuclei from the 
surface of the neutron star. This is caused by the lack of ef- 
ficient cascading in the inner pulsar magnetosphere which 
products impinge on the region of the polar cap and al- 
low efficient extraction of the iron nuclei from the neutron 
star surface. The cascades stop developing in the pulsar 
magnetosphere when the pulsar period becomes too low 
due to the pulsar rotational energy losses. This happens 
when pulsars reach the age of ~ 10 6 — 10 yrs, depending 
on the value of their surface magnetic field. Note however, 
that Giller & Lipski (2002) explain this low energy cut-off 
by the decay of the surface magnetic field of the pulsars 
on a typical time scale of ~ 10 7 yrs. 

The models (i) and (iii) for the initial parameters of 
the pulsar population are not consistent with the observed 
spectral shape of the CRs, since they give too narrow spec- 
tra with the peak somewhere between the knee and the 
ankle. Model (ii), postulating a correlation between the 
initial pulsar period and the strength of the surface mag- 
netic field of the the pulsar (Xu et al. 2001), overproduce 
CRs in the region of the ankle by about an order of mag- 
nitude, even for the assumed here relatively strong depen- 
dence of the lifetime of particles inside the Galaxy on their 
energy. 

The numerical calculations presented here allow us to 
include the energy loss processes of nuclei, accelerated 
by the pulsar, during their propagation inside the PWNe 
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and also consider their escape conditions from the PWNe. 
Therefore, we can predict the mass composition of nuclei 
injected from the PWNe into the Galaxy and compare it 
with the observations. The mass composition estimated in 
terms of our best fit model (iv) shows rise in the energy 
range ~ 10 6 — 10 7 GeV, and slower fall between ~ 10 8 GeV 
and a few 10 9 GeV, provided that the light extragalactic 
component expected at the highest energies extends down 
to 3 x 10 7 GeV. This behaviour of the mass composition is 
in general agreement with the most observational results 
available in the literature and collected in the paper by 
Horandcl (2004) (see Fig. 4). However, if the extragalac- 
tic component breaks at 10 9 GeV, then our model predicts 
heavy composition below 10 9 GeV as suggested by some 
other experiments (Dawson et al.1998, Ave et al. 2003, 
Dova et al. 2003). 

The models considered by us have also other interest- 
ing, worth to be mentioned features. For example, model 
(ii), postulating that a significant part of pulsars are born 
with short initial periods and strong surface magnetic 
fields (correlation between these two parameters) , predicts 
domination of the mass composition at energy ~ 10 9 GeV 
by protons from decay of neutrons released from nuclei 
inside the PWNa. Note that the existence of an addi- 
tional pulsar population inside the Galaxy with extreme 
parameters seems to be required by the observations of 
anomalous X-ray pulsars and magnetars. It can be re- 
sponsible for the observed anisotropics of CRs from direc- 
tions of the Galactic Centre and the Cygnus region (see 
Bednarek 2002, Bednarek 2003). This pulsar population 
might be also responsible for the highest energy compo- 
nent of the CR spectrum (> 10 19 eV), provided that heavy 
nuclei are able to escape from the PWNe (see Blasi et 
al. 2000, and Fig. Id in this paper). Such model predicts 
very heavy composition of the highest energy CRs. 

It might be surprising that we postulate a similar or- 
der contribution to the cosmic ray spectrum around the 
knee region from particles accelerated at the outer super- 
nova shock waves and from the pulsars. However, it is 
likely that the initial parameters of pulsars and expand- 
ing supernova envelopes, which originate in this same phe- 
nomenon, are in some way related. Pulsars with more ex- 
treme parameters (initial periods, surface magnetic fields) 
seem to be produced by explosions of type Ib/c super- 
novae which progenitors rotate fast and have relatively 
smaller envelopes just before explosion. Therefore, their 
envelopes have larger velocities, expanding to larger dis- 
tances and producing outer shocks with larger dimensions 
than in the case of more massive type II supernovae. Such 
shocks allow more efficient particle acceleration. Moreover, 
the exact shape of the cosmic ray spectrum in the knee 
region is not precisely known leaving place for some spec- 
ulations concerning its smoothness. For example, Erlykin 
& Wolfendale (1997) argue for the existence of fine struc- 
tures in the knee region of the cosmic ray spectrum to 
be interpreted as contribution from nuclei with different 
mass numbers. 



The rate of pulsar formation inside the Galaxy, es- 
timated on - 1/(20 - 250 yr) (Lyne et al. 1985; 
Lorimer 2003, Vranesevic et al. 2003), is about three or- 
ders of magnitude higher than the GRB event rate in 
our Galaxy. These pulsars can supply hadrons with the 
rate required by the observed density of the CRs in the 
Galaxy. Therefore, the population of classical radio pul- 
sars observed in the Galaxy can provide at least similar 
contribution to the galactic CRs as recent GRB model 
considered by Wick et al. (2004). 
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